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Abstract 

The proposed explanations are provided for the one-dimensional 
diffusion process with constant drift by using forward Fokker-Planck 
technique. We present the exact calculations and numerical evaluation 
to get the outflow probability in a finite interval, i. e. first passage time 
probability density distribution taking into account reflecting bound- 
ary on left hand side and absorbing border on right hand side. This 
quantity is calculated from balance equation which follows from con- 
servation of probability At first, the initial-boundary-value problem 
is solved analytically in terms of eigenfunction expansion which relates 
to Sturm-Liouville analysis. The results are obtained for all possible 
values of drift (positive, zero, negative). As application we get the 
cumulative breakdown probability which is used in theory of traffic 
flow. 

PACS numbers: 02.60.Lj, 02.50.-r, 02.50.Fz, 02.50. Ga 



1 Introduction 

Nowadays, the natural sciences deal with objects which have nondetermin- 
istic behaviour. Their descriptions can be found in theory of stochastic pro- 
cesses as a branch of probability theory There are different theoretical 
approaches for similar investigations using language of stochastic trajectories 




Figure 1: Schematic picture of the boundary-value problem showing the 
probability density p(x , t) in the interval a < x < b. 

as well as probability distributions. The fundamental equation which gives 
us the probabilistic description is the Fokker-Planck equation [3JIB]. Our mo- 
tivations for theoretical investigations in this field are given by application 
of the models of many-particles system which are considered in theoretical 
physics, i. e. physics of traffic flow [§]. Here we would like to find an ana- 
lytical solution for the special case when the stochastic variable belongs to a 
finite interval in terms of probability density distributions as well as cumu- 
lative probability [7j. The interval is defined as closed on left hand side and 
opened on right hand side. Due to these properties we introduce boundary 
conditions which determine the behaviour of the solution. Another impor- 
tant analysed quantity is the outflow (or breakdown) probability at right 
border which is found from the solution of Fokker-Planck equation by using 
balance equation. 

Let us consider the initial-boundary-value-problem (shown schematically 
in Fig. with constant diffusion coefficient D and constant drift coefficient 
v. Our task is to calculate the probability density p(x, t) to find the system 
in state x (exact in interval [x; x + dx\) at time moment t. The dynamics of 
p(x, t) is given by the forward drift-diffusion-equation as well as initial and 
boundary conditions [3 by following dynamics 
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dp(x, t) dp(x, t) -Q 9 2 p(x, t) 



dt dx 
^ dp(x,t) dj(x,t) 



dx 2 

-- 

dp(x, t) 
dx 



dt dx 
with flux j( x ,t) = vp(x,t) — D 

with the initial condition 

p(x, t = 0) = S(x — x ) , 
and two boundary conditions i. e. reflecting boundary at x = a 

dp(x, t) 



j(x — a,t) — vp(x = a,t) — D 

and absorbing boundary at x = b 

p(x — b, t) — . 



dx 







(1) 
(2) 
(3) 

(4) 
(5) 
(6) 



It is convenient to formulate the drift-diffusion problem in dimensionless 
variables. For this purpose we define new variables y and T by 



y 



x — a 
b — a 



and T 



D 



(6-o) : 



t 



As a result, the system of partial differential equations (fT|) 
rewritten as 

dP(y, T) _ n dP(y,T) + d 2 P(y, T) 



dT 



dy 



dy 1 



with initial condition 

P(y,T = 0) = 6(y-y ) 
reflecting boundary at y = 



J(y = 0,T)=QP(y = 0,T) 



dP(y,T) 



dy 



y=o 



and absorbing boundary at y = 1 

P(y = l,T)=0 



(7) 
can be 

(8) 
(9) 

(10) 
(11) 
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Hence, our problem has only one dimensionless control parameter £7 = 

v 

— (b — a) (scaled drift v which may have positive, zero, or negative val- 
ues). The parameter Q has the same meaning as Peclet number which has 
been used in j^j. 

The system of equations (JBJ) - (JUJ) will be solved exactly by applying the 
forward technique The main idea is to obtain the solution of Fokker- 
Planck equation and after that the first passage time distribution in terms 
of probability density. Both quantities will be presented as eigenfunction 
expantions. The survival probability and moments of first passage time can 
be calculated differently by using backward drift-diffusion equation. These 
results shown in [TJ |21 E] do not give the complete solution of the problem 
under consideration. Our presented analysis of the reference system (jHJ) - 
(fTTj) is the key result in order to study more complicated situations with 
nonlinear drift function 



2 Solution in terms of orthogonal eigenfunc- 
tions 

To find the solution of the well-defined drift-diffusion problem, first we take 
the dimensionless form (jHJ) - and use a transformation to a new function 
Qby 

Q(y,T) = e-^P(y,T). (12) 

This results in a dynamics without first derivative called reduced Fokker- 
Planck-equation 

^ = 4«^ + ^" ( 13 ) 
According to (fT2"|) the initial condition is transformed to 

Q(y,T = 0) = e-^°P(y,T = 0) , (14) 
whereas the reflecting boundary condition at y — becomes 



?g ( ^o,T) "tort 



2 ™ ' dy 



= 0, (15) 
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and the absorbing boundary condition at y = 1 now reads 

Q(y=l,T) = 0. (16) 

The solution of reduced equation (JT3j) can be found by the method of 
separation of variables p. Making a separation ansatz Q(y,T) = x(T)ip(y), 
we obtain 

1 d X (T) _ 1 dV(y) 

Both sides should be equal to a constant. This constant is denoted by —A, 
where A has the meaning of an eigenvalue. The eigenvalue A should be real 
and nonnegative. 

Integration of the left hand side gives exponential decay 

X (T)= Xo exp{-AT} (18) 

with x(T = 0) = Xo an d setting xo = 1- 

Let us now define the dimensionless wave number ka.sk 2 = A. The 
right-hand side of eq. ()17j) then transforms into the following wave equation 



dy 2 V 4 

Further on, we introduce a modified wave number k 2 = k 2 — Q 2 /4. Note that 
k = +\/k 2 — f2 2 /4 may be complex (either pure real or pure imaginary). 

First we consider the case where k is real. A suitable complex ansatz for 
the solution of the wave equation (fT§j) reads 

ij)(y) = C* exp{+iky} + C exp{-iky} (20) 

with complex coefficients C = A/2 + i B /2 and C* = A/2 — i B/2 chosen in 
such a way to ensure a real solution 

ip(y) = A cos(%) + B sin (fa/) . (21) 

The two boundary conditions (|T3|l and (fTtjj) can be used to determine the 
modified wave number k and the ratio A/B. The particular solutions are 
eigenfunctions i/j m (y), which form a complete set of orthogonal functions. As 
the third condition, we require that these eigenfunctions are normalised 



tinWv = 1 • (22) 
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In this case all three parameters k, A, and B are defined. 
The condition for the left boundary (|15|) reads 



%(y = o) - ^ 

2 ; dy 

After a substitution by (}2T?|) it reduces to 







or 



—A = kB . 
2 

The condition for the right boundary (fTo^) 

V(y = i) = o 



f23) 



(24) 
(25) 

(26) 



gives us 



or 



C* exp{-M&} + Cexp{-ifc} = 



A cos ( A; I + B sin ( A; 







(27) 



(28) 



By putting both equalities (f23J) and (|2*H|) together and looking for a non- 
trivial solution, we arrive at a transcendental equation 

%— (exp{+ik} — exp{— ik} j = k (exp{+ik} + exp{— ik} j (29) 



or 



sin ( k ) + k cos ( k 



respectively 



tan f k 



(30) 



(31) 



which gives the spectrum of values k m with m — 0, 1, 2, . . . (numbered in such 
a way that < k < k\ < k 2 < . . .) and the discrete eigenvalues X m > . 
Due to (|2*Tj) and (|2~Kj) . the eigenfunctions can be written as 



^m(y) = Rm cos ( ) sin ( k m ) - cos ( fc m ) sin ( k m y 



(32) 
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where R m = A m j sin ( k 



-B m / cos ( k, 



identity sin (a — (3) = sin a cos (3 — cos a sin f3, eq. 

i) m {y) = R m sm k m (l - y) 



Taking into account the 
reduces to 

(33) 



The normalisation constant R m is found by inserting (|33|) into (J22)) . Calcula- 
tion of the normalisation integral by using the transcendental equation (j3*Uj) 
gives us 



R, 



sin 



R 2 , 



k m (1 - y) 



n 



2 ki + n*/4 



= 1 , 



sin ( 2k 



(34) 



and hence (j3l?|) becomes 



1 + 



n _ i 

2 ^+^2/4 



■ sm 



k m (l - y) 



or 




sm 



- n 2 /4 (i - y) 



(35) 



(36) 



This calculation refers to the case Q > —2 where all wave numbers k m or 
k m = y/k^ — VL 2 /A are real and positive. 

However the smallest or ground-state wave vector k vanishes when Q 
tends to —2 from above, and no continuation of this solution exists on the 
real axis for Q < —2. A purely imaginary solution k = in appears instead, 
where kq is real, see Fig. El In this case (for Q < —2) a real ground-state 
eigenfunction ipo(y) can be found in the form (|2()jl where C = A/2 + B/2 and 
C* = A/2 - B/2, i. e., 



ipo{y) = Acosh(K y) + B sinh(/c Z/) ■ 



(37) 



The transcendental equation for the wave number k = ik, can be written 
as the following equation for k 



■ sinh (kq) + kq cosh (kq) = . 



(38) 
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Figure 2: The wave number k (Q > — 2) respectively fto (^ < — 2) and 
eigenvalue Ao for ground state m — 0. The thin straight line shows the 
approximation k,q « — Q/2 valid for large negative Q < — 5. 

As compared to the previous case Q > — 2, trigonometric functions are 
replaced by the corresponding hyperbolic ones. Similar calculations as before 
yield 



Note that k,q = —iko is the imaginary part of ko and = — k$. As regards 
other solutions of (|3T)|) called excited states, i. e., those for k m with m > 0, 
nothing special happens at Q = —2, so that these wave numbers are always 
real. The situation for ground state m = at different values of dimensionless 
drift parameter Q is summerized in Tabled which presents the solutions hq 
from transcendental equation (pIHj) together with A = — + ^ 2 /4 and k 
from transcendental equation (JHUj) together with eigenvalues A = k$ +f2 2 /4. 
Table|21shows the behaviour of lowest wave numbers k m with m = 0,1, ... ,5. 
The results are plotted in Fig. El 

In general (for arbitrary Q), the eigenfunctions are orthogonal and nor- 
malised, i. e., 



Jo 

Figure 0] shows the ground eigenstate (m = 0) for different parameter values 
Q, whereas Fig. Ogives a collection of eigenstate functions (m — 0, 1, . . . , 5) 




(39) 




(40) 



for n = -5.0 and = 3.0. 
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Table 1: The ground-state wave number kq (for Q < —2) and ko (for Q > —2) 
and eigenvalue Ao depending on the dimensionless drift parameter Q. 





K 


Ao 


— q nn 


A AQQ 


n m n 


—R ^n 


A 9/L8 


U.U 1 o 


—8 nn 


3 QQ7 


n n2i 


-7.50 


3.745 


n.n3i 


-7.00 


3.493 


n.n45 


-6.50 


3.240 


n.n64 


-6.00 


2.984 


0.091 


-5.50 


2.726 


0.128 


-5.00 


2.464 


0.178 


-4.50 


2.195 


0.245 


-4.00 


1.915 


0.333 


-3.50 


1.617 


0.446 


-3.00 


1.288 


0.591 


-2.50 


0.888 


0.774 


-2.00 


0.000 


1.000 





k 


Ao 


—9 nn 


n nnn 


1 nnn 


—i ^n 


n 84^ 




—1 nn 


1 1 65 


1 608 


-0.50 


1.393 


2.004 


0.00 


1.571 


2.468 


0.50 


1.715 


3.005 


1.00 


1.836 


3.623 


1.50 


1.939 


4.325 


2.00 


2.928 


5.116 


2.50 


2.106 


5.999 


3.00 


2.174 


6.979 


3.50 


2.235 


8.058 


4.00 


2.288 


9.239 


4.50 


2.337 


10.525 


5.00 


2.381 


11.917 



Table 2: The wave numbers k m (m = 0, 1, . . . , 5) depending on the dimen- 
sionless drift parameter Q. 





-10.0 


-5.0 


-2.0 


-1.0 


0.0 


1.0 


2.0 


5.0 


10.0 


m = 


4.999 


2.464 


0.000 


1.165 


1.571 


1.836 


2.028 


2.381 


2.653 


m = 1 


3.790 


4.172 


4.493 


4.604 


4.712 


4.816 


4.913 


5.163 


5.454 


m = 2 


7.250 


7.533 


7.725 


7.789 


7.854 


7.917 


7.979 


8.151 


8.391 


m = 3 


10.553 


10.767 


10.904 


10.949 


10.995 


11.040 


11.085 


11.214 


11.408 


m = 4 


13.789 


13.959 


14.066 


14.101 


14.137 


14.172 


14.207 


14.310 


14.469 


m = 5 


16.992 


17.133 


17.220 


17.249 


17.279 


17.308 


17.336 


17.421 


17.556 
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Figure 3: The parameter dependence of wave numbers k m (Q) and eigenvalues 
A m (f2) for ground state m = and excited states m — 1, 2, 3. 
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Figure 5: The eigenfunctions if) m (y) for m — 0, 1, 2, 3 and for = —5.0 (left) 
and VL = 3.0 (right). 

In the following, explicit formulae (where ip m (y) is specified) are written 
for the case Q > — 2. 

In order to construct the time-dependent solution for Q(y,t), which ful- 
fills the initial condition, we consider the superposition of all particular so- 
lutions with different eigenvalues A m 

oo 

Q{y,T) = Y,C m e- XmT ^ m {y). (41) 

m=0 

By inserting the initial condition 

P(y, T = 0) = e^Q(y, T = 0) = 5(y - y ) (42) 
into (}4*Tj) we obtain 

oo 

Y,C m ^m(y) = e-% y 5(y-y ) . (43) 

m=0 

Now we expand the right hand side of this equation by using the basis of 
orthonormalised eigenfunctions (|35|) and identify C m with the corresponding 
coefficient at if) m , i. e., 

C m = e% y 5(y - y^Andy = e-% yo il> m (y ) • (44) 
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This allows us to write the solution for P(y,T) as 



P{y,T) = e 



K2/-2/0) 



m=0 



- A m T 



?Pm(yo)?Prn(y) 



with eigenfunctions (J35J) and (J39"j) of ground state (m = 0) 

feo(l-y) 



^o(y) 



1 + 



n 1 



sin 



2 

>/3 (1-2/) , 



ft > -2 



n = -2 



sinh [« (1 — j/)] , ft < - 2 



2 - K g+n 2 /4 
and all other eigenfunctions (p?B*j) 



1 + 2 



sin 



fc m (l -y) 



TO = 1,2, 



2 fc^+na/4 

The eigenvalue of ground state (to = 0) is given by 

' ~k 2 + ft 2 /4 , ft > - 2 



1, ft = -2 

k - k§ + ft 2 /4 , ft < - 2 



and all others are 



A, 



(45) 



(46) 



(47) 



(48) 



~k 2 m + n 2 /4 m = l,2,..., (49) 

where the wave numbers are calculated from transcendental equation (pTT|) 

2 ~ 

k : tan k = — — k ft > — 2 (50) 

(51) 
(52) 



tanh k = — ^ K o ft < — 2 



tan k r . 



TO = 1, 2 



The set of Figures El illustrates the time evolution of probability density (J45|) 
choosing different parameter values ft. 



12 





2- 



H 



T = 0.01 
T = 0.05 
T = 0.1 
T = 0.5 



1„ 



T = 0.01 
• T = 0.05 
T = 0.1 
T = 0.5 



0.5 

y 



0.5 

y 



Figure 6: The solution of drift-diffusion Fokker-Planck equation with initial 
condition yo = 0.5 for different values of the control parameter Q, i. e. Q = 
-5.0 (top left), Q — -2.5 (top right), ft = 0.1 (bottom left), = 3.0 
(bottom right). 
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3 First passage time probability density 



It has been shown in previous sections that the probability density P(y, T) is 
not normalized under given restrictions, i. e. reflected at y = and absorbed 
at y = 1. Due to that fact, let us apply here the balance equation in the 
open system given in dimensionless variables by 



V(T,y 



d_ 



P(y,T) dy 



(53) 



which relates the probability P(y, T) that the system is still in a state 
y G [0,1] with the probability flux V(T,y = 1) out of this interval at the 
right absorbing boundary y = 1 at time moment T. Hence, V(T,y = 1) is 
the first passage time probability density [H 12 El • It can be calculated by 
using obtained results of previous section. The first passage time probabil- 
ity density distribution V (breakdown probability density) depending on VL 
reads as follows 



I. VL > -2 

CO 

V(T,y = l) = 2e^^J2j 



m=0 



+ 



n 



sin 



2. n = -2 



(i-i/d) 



V(T,y = l) = e 

CO 

Ee 
7 



2 kl+ny* 
3(l-y )e- T 

(fc^+l)T_ 

km Sill 



2/o J 



(54) 



1 ,v m 



m=l " 

3. n < -2 

V(T,y = l) = 2e§ (1 " w) 

e -(- K 2+C 2 /4)T 



fc m (l - y ) 



(55) 



x 



1 + 



2 - K g+Q 2 / 4 



K sinh [«o(l - y )] 



+ E i | n i — fcmsin 

TO =1 1 "r 2 



- 2/0 ) 



(56) 



14 



Figure 7: The first passage time probability density distribution V(T, y — 1) 
for Vt < -2 (left) and Vt > -2 (right). 

The outflow distribution V(T, y — 1) is shown in Fig. [7| (with different values 
of dimensionless drift f2) as well as in Fig. |H] (with different values of initial 
condition y ). 

4 Cumulative breakdown probability 

The probability that the absorbing boundary y — 1 is reached within certain 
observation time interval < T < T 6 S is given by the cumulative (break- 
down) probability 



with V(T,y = 1) from For T f, s — > oo we have W ^ 1. Generally, we 

obtain 

1. > -2 




(57) 



o 



W(n,T obs ) = 2 e §( 1 -^^ 



1 _ e -(^+^ 2 /4)T o6s 



/c m sin fc m (l-y ) 



m=0 



+ fi 2 /4 + n/2 



(58) 
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Figure 8: Short time behaviour of first passage time probability density 
distribution V(T, y = 1) for different initial conditions < yo < 1 showing 
time lag. 



2. n 



(i-yo) 



W(n,T obs ) = e 



3. VL < -2 

w(n,T obs ) 



3 (1 - e~ T °*°) (1 - y ) 
1 - e 



m=l 



sin 



fem(l - 2/0 ) 



• (59) 



2e2 



"(1-2/0) 



1 - e 



^ + n 2 /4 + ^/2 



K sinh [«o(l - y )] (60) 



~ 1 _ e -(kl+n 2 /i)T obs _ 

+ > -= fc m sin 

£l kl + W/4 + n/2 



k m {l - y ) 



Figure El shows W (fl,T b s ) as a function of observation time T t, s (left) as 
well as parameter dependence Q (right). 
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Figure 9: The probability W(Q,T b s ) ()57|) as function of observation time 
T bs with fixed fl (left) and vice versa (right). 

5 Limit case for large positive values of the 
control parameter 

Consider parameter limit Q — > +00 which corresponds either to large positive 
drift v and/or large interval b — a, or to a small diffusion coefficient D. In 
this case, for a given m, the solution of the transcendental equation can be 
found in the form k m = 7i(m + 1) — e m , where e m is small and positive. From 
the periodicity property we obtain 

cosfc m = cos(7r(m + 1) - e m ) = -(-l) m cos{e m ) = -(-l) m + O (e 2 m ) 

sink m = sin(7r(m + 1) - e m ) = (-l) m sin(e m ) = (-l) m e m + O (e 3 m ) . 
By inserting this into the transcendental equation (JHUj) . we obtain 

e m = ^n(m + l) + 0(n- 2 ) , (61) 
sm(k m ) = ^(-l) m vr(m+l) + 0(fi- 2 ) . (62) 
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In this approximation the normalisation integral for large f2 and the initial 
condition y — > can be written as 

7_ , ™ k m sin ( k m ) 

I = jv(T,y = l) d T = 2e^^-^ (63) 

n/2 -4 (-l)-(vrm) 2 _ /2 ~ -2 (-l^m) 2 
^n7r 2 m 2 + fi 2 /4 ^ fivr 2 m 2 + ft 2 /4' 

m=l ' m=— oo ' 

Further on we set (— l) m = e t7TTn and, in a continuum approximation, replace 
the sum by the integral 

oo 

0/9 f -2 e i7rm (7rm) 2 , , s 

J Q Tc 2 m 2 + 9? /A v ' 

— oo 

Now we make an integration contour in the complex plane, closing it in 
the upper plane (Imm > 0) at infinity where \e lnm \ is exponentially small. 
According to the residue theorem, it yields 

/ = 2ni Res(mj) = 27riRes(m ) , (65) 

i 

where mo = lp is the location of the pole in the upper plane, found as a 
root of the equation 7r 2 m 2 + f2 2 /4 = 0. According to the well-known rule, 
the residue is calculated by setting m = mo in the enumerator of (JMj) and 
replacing the denominator with its derivative at m — tuq. It gives the desired 
result I — 1, i. e., the considered approximation gives correct normalisation 
of outflow probability density V(T,y = 1) at the right boundary. 

The probability distribution function P(y,T) given by (j43|) can also be 
calculated in such a continuum approximation. In this case the increment of 
wave numbers is 

Ak rn = k m+1 -k m = n + e m - e m+1 ~ w (l - ~ - ^/q, • (66) 

Note that in this approximation for Q —>■ oo the normalisation constant R m 
in f!34|) is related to the increment Ak via 

^ = fTi?^= = TT2jn H l A ' k - ■ (67) 
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Hence, the equation (J45j) for the probability density can be written as 

CO 

P(y, T) = 2e^°) ^ R 2 m e~ x - T sin \k m (1 - y )] sin \k m (1 - y)] (68) 



m=0 



7T 



e f dr-Kd) e-(' k - +n2 ^ T sin [fc m (1 - y )] sin [fc m (1 - y) 



m=0 



In the continuum approximation we replace the sum by the integral 



P(y, T) ~ l e T(»-w) / e -(k 2 +^/A)T gin k (1 _ y(j )l gin \~k (1 - y)] dk (69) 

7T 



7T 



* (y - yo) 



cos 



k(2-y-yo) 



dk . 



In the latter transformation we have used the identity 
sin a sin (3 = | (cos(a — /3) — cos(a + (3)). The resulting known integrals 
yield 

1 n. 



P(y,T) 



VAttT 



(y-vo) (z-y-yp) 

g 4T — g 4T 



(70) 



The approximation (J7U|) is shown in Fig. El For short enough times 4T <C 
(2 — y — y ) 2 the second term is very small. Neglecting this term, Eq. (J7(Jj) 
reduces to the known exact solution for natural boundary conditions. 
Based on (J7UJ), it is easy to calculate the probability flux 



J{y,T) = nP{y,T)--?-P(y,T) 



(71) 



and the first passage time distribution V(T) = J(y = 1, T) which takes a 
particularly simple form 



V(T) 



1 — W (i-yp-OT) 2 

-e 4T 



The cumulative breakdown probability ()57|) is then 

T 

± obs 

W(Q,T = T obs )= J -j==e - dT. 



(72) 



(73) 



19 




Figure 10: Comparison of probability density P(y,T) in drift-diffusion- 
dynamics with finite boundaries for two time moments. Parameter value is 
Q = 3.0; initial condition is yo = 0.5. The solid lines represent the exact 
result (jlSJ); dotted lines display the approximation (J7UJ). 



6 Relationsship to Sturm— Liouville theory 



The particular drift-diffusion-problem over a finite interval with reflecting 
(left) and absorbing (right) boundaries belongs to the following general math- 
ematical theory named after Jacques Charles Francois Sturm (1803-1855) 
and Joseph Liouville (1809-1882). 

The classical Sturm-Liouville theory considers a real second-order linear 
differential equation of the form 



d 

dx 



' dip' 
P{x) di- 



+ q(x)ip = Xw(x)ip 



(74) 



together with boundary conditions at the ends of interval [a, b] given by 

dip 



a\ip{x = a) + ct2 
(3 l iP{x = b)+(3 2 



dx 
dip 
dx 







. 



(75) 
(76) 



x=b 



The particular functions p(x),q(x),w(x) are real and continuous on the 
finite interval [a, b] together with specified values at the boundaries. The aim 
of the Sturm-Liouville problem is to find the values of A (called eigenvalues 
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X n ) for which there exist non-trivial solutions of the differential equation (|74j) 
satisfying the boundary conditions (|75j) and (J75J) . The corresponding solu- 
tions (for such A n ) are called eigenfunctions ip n {x) of the problem. 

Defining the Sturm-Liouville differential operator over the unit interval 
[0, 1] by 

d 



dx 



(77) 



and putting the weight it? (a;) to unity (it? = 1) the general equation (f71)l can 
precisely be written as eigenvalue problem 

Lil) = \ip (78) 

with boundary conditions (|75|)(a = 0) and (J75j) (6 = 1) written as 

B ip = Blip = . (79) 

Assuming a differentiable positive function p(x) > the Sturm-Liouville 
operator is called regular and it is self-adjoint to fulfil 



4>i ■ ^2 ■ 



(80) 



Any self-adjoint operator has real nonnegative eigenvalues Ao < Ai < . . . < 
A n < . . . — ► oo. The corresponding eigenfunctions ip n (x) have exact n zeros 
in (0, 1) and form an orthogonal set 



^ n (x)ip m (x)dx = 5 n 



[81) 



The eigenvalues A ra of the classical Sturm-Liouville problem (|73j) with 
positive function p(x) > as well as positive weight function w(x) > 
together with separated boundary conditions (|75j) and (fTUj) can be calculated 
by the following expression 



A n / ll> n (x) 2 w(x)dx 



\p(x) (dip n (x) I dx) 2 + q(x)ipn{x) 2 ] dx 

b 

p(x)ip n (x) (dip n (x) I 'dx) . (82) 
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The eigenf unctions are mutually orthogonal (m ^ n) and usually normal- 
ized (m = n) 

ip n {x)ip m {x)w(x)dx = 5 mn (83) 

known as orthogonality relation (similar to ()81|). 

Comming back to the original drift-diffusion problem written in dimen- 
sionless variables over unit interval < y < 1 and recalling (|17l the separation 
constant A appears in the following differential equation 



(84) 



which can be related to the regular Sturm-Liouville eigenvalue problem via 
p(y) = 1 > 0; w(y) = 1 > and g(y) = fi 2 /4. 

The boundary conditions given by (J2*H|) and (J2B|) can be expressed as 



0.^ = o ) + ( -i).f 



1 ■ V(|/ 



0- 



y=0 



j/=l 









(85) 
(86) 



in agreement with (J75|) and (J75|) . 

The up-to-now unknown separation constant A has a spectrum of real 
positive eigenvalues which can be calculated using (}82j) from 



A, 



dif> n (y) 
dy 



fi 5 



dx — 



dip n {y) 



dy 



(87) 



taking into account normalized orthogonal eigenfunction (|83|) 



ip n (y)ipm(y)dy 



(88) 



7 Conclusions 

The presented paper shows the analytical method how to solve drift-diffusion 
initial-boundary-value problem for the case of reflecting and absorbing bound- 
aries [H 12 IH El EI • On the basis of Sturm-Liouville theory, the set of eigen- 
values with corresponding eigenfunctions has been found. Here we have paid 
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our attention to wave number calculations from transcendental equations. 
The equations have been solved numerically by Newton method. The main 
problem which has been solved was the dependence of obtained results on 
drift value, i. e. different cases of control parameter Q < —2, Q = — 2 and 
f2 > —2. First case of £7 < — 2 corresponds to the situation when it is difficult 
and probably impossible, with significant small probability and for long times 
only, to leave the interval due to the large negative value of drift. The case of 
Q = —2 has been considered as limit case and the corresponding solution has 
been found. The opposite case of Q > — 2 shows the usual situation when 
the system reaches the right border relatively fast. As application, the first 
passage time distribution as well as the cumulative probability have been 
calculated. The case of large positive values of fl has been investigated in 
detail and has been obtained as approximation. 
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